/**
 *
 *      @file PoseScore.c
 *      @brief Routines for the PoseScore and RankScore statistical potentials
 *
 *      @author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *      @date 01/02/2012
 *
 *      This program is free software; you can redistribute it and/or modify
 *      it under the terms of the GNU General Public License as published by
 *      the Free Software Foundation version 2 of the License.
 *
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU General Public License for more details.
 *
 */

#include "PoseScore.h"

float get_PoseScore_pair(int lig, int prot, float dist);
float get_RankScore_pair(int lig, int prot, float dist);
float get_RankScore_pair_ipol(int lig, int prot, float dist);
int translate_gaff_to_PoseScore(int gaff_type);


/**
 *
 *	@brief Get RankScore for a molecule
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *	@param lig MOL2 with the ligand
 *	@param prot MOL2 with the protein
 *	@return Score
 *
 */
float get_RankScore(MOL2 *lig, MOL2 *prot)
{
	int i = 0, j = 0, k = 0;
	float dx = 0.0f, dy = 0.0f, dz = 0.0f;
	float dist = 0.0f;
	float energy = 0.0f;

	for ( i = 0; i < lig->n_atoms; ++i) {
		for ( j = 0; j < prot->n_atoms; ++j) {
			dx = lig->x[i] - prot->x[j];
			dy = lig->y[i] - prot->y[j];
			dz = lig->z[i] - prot->z[j];
			dist = dx * dx + dy * dy + dz * dz;

			if ( dist < 36) /* 6A cutoff */
				energy += get_RankScore_pair_ipol( translate_gaff_to_PoseScore(lig->gaff_types[i]), prot->gaff_types[j], sqrt(dist));
		}
	}

	return energy;
}


/**
 *
 *	@brief Get PoseScore for a molecule
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *	@param lig MOL2 with the ligand
 *	@param prot MOL2 with the protein
 *	@return Score
 */
float get_PoseScore(MOL2 *lig, MOL2 *prot)
{
	int i = 0, j = 0, k = 0;
	float dx = 0.0f, dy = 0.0f, dz = 0.0f;
	float dist = 0.0f;
	float energy = 0.0f;

	for ( i = 0; i < lig->n_atoms; ++i) {
		for ( j = 0; j < prot->n_atoms; ++j) {
			dx = lig->x[i] - prot->x[j];
			dy = lig->y[i] - prot->y[j];
			dz = lig->z[i] - prot->z[j];
			dist = dx * dx + dy * dy + dz * dz;

			if ( dist < 36) /* 6A cutoff */
				energy += get_PoseScore_pair( translate_gaff_to_PoseScore(lig->gaff_types[i]), prot->gaff_types[j], sqrt(dist));
		}
	}

	return energy;
}

/**
 *
 *	@brief Rough translation between GAFF atom types and Pose/RankScore types
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *	@param gaff_type GAFF atom type
 *	@return PoseScore atom type
 *
 */
int translate_gaff_to_PoseScore(int gaff_type)
{
	int ret = -1;

	/* Carbon block */
	if ( gaff_type == C1)
		ret = 0;
	else if ( gaff_type == C2 || (gaff_type >= 23 && gaff_type <= 35) )
		ret = 1;
	else if ( gaff_type == C3 )
		ret = 2;
	else if ( gaff_type == CA )
		ret = 3;
	else if ( gaff_type == N1)
		ret = 5;
	else if ( gaff_type == N2 || (gaff_type >= 42 && gaff_type <= 46) )
		ret = 6;
	else if ( gaff_type == N3)
		ret = 7;
	else if ( gaff_type == N4)
		ret = 8;
	else if ( gaff_type == NA)
		ret = 9;
	else if ( gaff_type == N)
		ret = 10;
	else if ( gaff_type == NO)
		ret = 11;
	else if ( gaff_type == O)
		ret = 12;
	else if ( gaff_type == OH || gaff_type == OS || gaff_type == OW)
		ret = 13;
	else if ( gaff_type == S || gaff_type == S2)
		ret = 16;
	else if ( gaff_type == S6 || gaff_type == SH || gaff_type == SS)
		ret = 17;
	else if ( gaff_type == SX || gaff_type == SY)
		ret = 18;
	else if ( gaff_type == S4)
		ret = 19;
	else if ( gaff_type <= 66 && gaff_type >= 57)
		ret = 21;
	else if ( gaff_type == F)
		ret = 22;
	else if ( gaff_type == CL)
		ret = 23;
	else if ( gaff_type == BR)
		ret = 24;
	else if ( gaff_type == I)
		ret = 25;

	return ret;
}


/**
 *
 *	@brief Interpolated PoseScore for a pair of atoms
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *	@param lig Ligand atom type
 *	@param prot Protein atom type
 *	@param dist Distance between atoms in A
 *	@return Score
 *
 */
float get_PoseScore_pair_ipol(int lig, int prot, float dist)
{

	int bin = 0,  bin2 = 0;
	float dist1 = 0.0f, dist2 = 0.0f;

	if ( dist > 5.9 || dist <= 0.1 || lig >= 26 || lig < 0 || prot >= 158 || prot < 0)
		return 0.0f;

	bin = (int)(dist / 0.1f);
	bin2 = bin + 1;

	dist1 = bin * 0.1f;
	dist2 = bin2 * 0.1f;

	return (((dist - dist1) / (dist2 - dist1)) * ( PoseScore_pot[lig][prot][bin2] - PoseScore_pot[lig][prot][bin])) + PoseScore_pot[lig][prot][bin];

}

/*
 *	@brief Get PoseScore for a pair of atoms
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *      @param lig Ligand atom type
 *      @param prot Protein atom type
 *      @param dist Distance between atoms in A
 *      @return Score
 */
float get_PoseScore_pair(int lig, int prot, float dist)
{

	int bin = 0;

	if ( dist > 6.0 || dist <= 0 || lig >= 26 || lig < 0 || prot >= 158 || prot < 0)
		return 0.0f;

	bin = (int)(dist / 0.1f);
	return PoseScore_pot[lig][prot][bin];
}


/*
 *     @brief Get RankScore for a pair of atoms
 *     @author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *     @param lig Ligand atom type
 *     @param prot Protein atom type
 *     @param dist Distance between atoms in A
 *     @return Score
 */
float get_RankScore_pair(int lig, int prot, float dist)
{

	int bin = 0;

	if ( dist > 6.0 || dist <= 0 || lig >= 26 || lig < 0 || prot >= 158 || prot < 0)
		return 0.0f;

	bin = (int)(dist / 0.1f);
	return RankScore_pot[lig][prot][bin];
}


/**
 *
 *      @brief Interpolated RankScore for a pair of atoms
 *      @author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *      @param lig Ligand atom type
 *      @param prot Protein atom type
 *      @param dist Distance between atoms in A
 *      @return Score
 *
 */
float get_RankScore_pair_ipol(int lig, int prot, float dist)
{

        int bin = 0,  bin2 = 0;
        float dist1 = 0.0f, dist2 = 0.0f;

        if ( dist > 5.9 || dist <= 0.1 || lig >= 26 || lig < 0 || prot >= 158 || prot < 0)
                return 0.0f;

        bin = (int)(dist / 0.1f);
        bin2 = bin + 1;

        dist1 = bin * 0.1f;
        dist2 = bin2 * 0.1f;

        return (((dist - dist1) / (dist2 - dist1)) * ( RankScore_pot[lig][prot][bin2] - RankScore_pot[lig][prot][bin])) + RankScore_pot[lig][prot][bin];

}




/**
 *	@brief Pose/Rank Score protein typing
 *	@author Alvaro Cortes Cabrera <acortes@cbm.uam.es>
 *	@param myres Name of the residue
 *	@param myname Name of the atom
 *	@return Atom type
 */
int get_PoseScore_protein_type(char *myres, char *myname)
{
	char *myname2 = NULL;
	int i = 0, j = 0, k = 0;
	int score = 0;
	int PoseScore_type = -1;

	myname2 = (char*)calloc(sizeof(char *), 3);
	j = 0;

	for ( i = 0; i < strlen(myname); i++) {
		if ( myname[i] >= 65 && myname[i] <= 122) {
			myname2[j] = myname[i];
			j++;
		}
	}

	myname2[j] = '\0';


	for ( j = 0; j < 158; j++) {
		if ( strstr(myres, PoseRank_prot_types_names[j][0]) != NULL) {
			score = 0;
			for ( i = 0; i < strlen(PoseRank_prot_types_names[j][1]); ++i)
				if ( PoseRank_prot_types_names[j][1][i] == myname2[i])
					++score;

			if ( score == strlen(PoseRank_prot_types_names[j][1]) && score == strlen( myname2))
				PoseScore_type = PoseRank_prot_types[j];

		}
	}

	free(myname2);

	return PoseScore_type;
}



